#*************         NHANES         *************
#*************      Analyze Code      *************
#*************    Cox 回归模型讲解    *************
#*************                        *************
#### 0.准备好环境 ####
library(tidyverse)
library(gtsummary)
library(tidyr) # drop_na 函数，用于快速去掉 NA
library(survey)
library(haven)
library(survival) # 生存函数
library(survminer) # 生存函数
library(readr) # 读取 mort 数据
library(dplyr) # 读取 mort 数据
library(plyr) # 读取 mort 数据
# 运行一下读取 mort 数据的函数
setwd("G:/BaiduNetdiskDownload/NHANES")
NHSMortRead <- function(file.path){
  srvyin <- file.path   # full .DAT name here
  srvyout <- "data" # shorthand dataset name here
  
  # Example syntax:
  #srvyin <- paste("NHANES_1999_2000_MORT_2019_PUBLIC.dat")   
  #srvyout <- "NHANES_1999_2000"      
  
  
  # read in the fixed-width format ASCII file
  dsn <- read_fwf(file=srvyin,
                  col_types = "iiiiiiii",
                  fwf_cols(seqn = c(1,6),
                           eligstat = c(15,15),
                           mortstat = c(16,16),
                           ucod_leading = c(17,19),
                           diabetes = c(20,20),
                           hyperten = c(21,21),
                           permth_int = c(43,45),
                           permth_exm = c(46,48)
                  ),
                  na = c("", ".")
  )
  
  # NOTE:   SEQN is the unique ID for NHANES.
  
  # Structure and contents of data
  str(dsn)
  
  
  # Variable frequencies
  
  #ELIGSTAT: Eligibility Status for Mortality Follow-up
  table(dsn$eligstat)
  #1 = "Eligible"
  #2 = "Under age 18, not available for public release"
  #3 = "Ineligible"
  
  #MORTSTAT: Final Mortality Status
  table(dsn$mortstat, useNA="ifany")
  # 0 = Assumed alive
  # 1 = Assumed deceased
  # <NA> = Ineligible or under age 18
  
  #UCOD_LEADING: Underlying Cause of Death: Recode
  table(dsn$ucod_leading, useNA="ifany")
  # 1 = Diseases of heart (I00-I09, I11, I13, I20-I51)
  # 2 = Malignant neoplasms (C00-C97)
  # 3 = Chronic lower respiratory diseases (J40-J47)
  # 4 = Accidents (unintentional injuries) (V01-X59, Y85-Y86)
  # 5 = Cerebrovascular diseases (I60-I69)
  # 6 = Alzheimer's disease (G30)
  # 7 = Diabetes mellitus (E10-E14)
  # 8 = Influenza and pneumonia (J09-J18)
  # 9 = Nephritis, nephrotic syndrome and nephrosis (N00-N07, N17-N19, N25-N27)
  # 10 = All other causes (residual)
  # <NA> = Ineligible, under age 18, assumed alive, or no cause of death data available
  
  #DIABETES: Diabetes Flag from Multiple Cause of Death (MCOD)
  table(dsn$diabetes, useNA="ifany")
  # 0 = No - Condition not listed as a multiple cause of death
  # 1 = Yes - Condition listed as a multiple cause of death
  # <NA> = Assumed alive, under age 18, ineligible for mortality follow-up, or MCOD not available
  
  #HYPERTEN: Hypertension Flag from Multiple Cause of Death (MCOD)
  table(dsn$hyperten, useNA="ifany")
  # 0 = No - Condition not listed as a multiple cause of death
  # 1 = Yes - Condition listed as a multiple cause of death
  # <NA> = Assumed alive, under age 18, ineligible for mortality follow-up, or MCOD not available
  
  # Re-name the dataset, DSN, to the short survey name then remove other R objects
  assign(paste0(srvyout), dsn)
  rm(dsn, srvyin, srvyout)
  # View(data)
  return(data)
}


#### 1. 常规Cox模型 ####
##### 1.1 KM & Logrank #####
head(lung)
# time-时间: 以天为单位的生存时间
# status-状态: 1 =删失，2 = 死亡

# survfit 用来计算 Kaplan-Meier 生存估计，它的主要论点包括:
# 使用函数 Surv ()创建的生存对象，数据集

fit <- survfit(Surv(time, status) ~ sex, data = lung)
fit
# fit 显示生存曲线的简短摘要。包括观察数、事件数、中位生存期和中位置信限区间。
# 具体的事件发生的每个时间点的信息，risk 数量、事件数量、删失数量、生存概率、生存概率的 95% 置信区间
# 可以提取每个时间点其中的特定信息
output <- data.frame(time = fit$time,
                     n.risk = fit$n.risk,
                     n.event = fit$n.event,
                     n.censor = fit$n.censor,
                     surv = fit$surv,
                     upper = fit$upper,
                     lower = fit$lower)
# 可视化结果
ggsurvplot(fit,
           conf.int = TRUE,
           risk.table.col = "strata", # Change risk table color by groups
           ggtheme = theme_bw(), # Change ggplot2 theme
           palette = c("#E7B800", "#2E9FDF"),
           xlim = c(0, 600))

# 比较生存曲线的对数秩检验-图中的 p value 的含义
surv_diff <- survdiff(Surv(time, status) ~ sex, data = lung)
surv_diff

# 使用多种因素的组合来计算 KM 生存曲线，但是本质上还是单个因素
# 只是取值是多个变量的取值合并
fit2 <- survfit(Surv(time, status) ~ sex + rx + adhere,
                data = colon)

ggsurv <- ggsurvplot(fit2, fun = "event", conf.int = TRUE,
                     ggtheme = theme_bw())

ggsurv$plot + theme_bw() + 
  theme (legend.position = "right")+
  facet_grid(rx ~ adhere)


##### 1.2 Cox 比例风险模型 #####
# Kaplan-Meier 曲线和 Logrank 检验,是单变量分析的例子。他们根据调查中的一个因素来描述生存情况，却忽略了其他因素的影响。
# 但是我们非常建议大家从单变量入手，结合可视化，对数据建立初步、扎实的理解

res.cox <- coxph(Surv(time, status) ~ sex, data = lung)
res.cox
summary(res.cox)
# 解读1：变量是否有统计学意义，“ z”的列给出 Wald 统计值。
# Wald 统计量评价给定变量的 β-coef 是否等于0，即 HR-exp(coef)是否等于1
# 它对应于每个回归系数与其标准误差(z = coef/se (coef))的比值。

# 解读2: 变量回归系数的符号 & 值
# Cox 模型结果中需要注意的第二个特征是回归系数(coef)的符号。
# 女性为2，男性为1，Cox 模型的汇总给出了第二组相对于第一组的风险比(HR) ，即女性对男性。
# 在这些数据中，性别的 β 系数 = -0.53, HR <1 表明女性的死亡风险(存活率)低于男性。
lung$sex #当性别变量变“大”的时候，即从男性到女性，风险降低，假设男性风险是 1，女性是 0.59 倍，降低了 41%

# 解读3: 模型的 p value
# 模型的全局统计显著性有3个不同的检验，当数据量较大的时候，可以任选，都差不多
# 数据量较小，建议选择似然比的结果- Likelihood ratio test，针对样本会更稳定

# 多元 Cox 也是一样的函数
res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data =  lung)
summary(res.cox)

##### 1.3 使用 gtsummary 优化输出 #####

m1_tbl_1 <-
  tbl_regression(
    res.cox
  )
m1_tbl_1

#### 2. 加权 Cox 模型 ####
##### 2.1 提取数据模块 ##### 
demo.i <- read_xpt("2015-2016/Demographics/demo_i.xpt")#要提前设置好数据存储的路径
colnames(demo.i)
dim(demo.i) # 9971

# 提取生存数据 

# 连接 National Death Index (NDI) ，补充 NHANES 参与人员的死亡数据
# 得到公开的 linked mortality files (LMF) 提供从调查参与日期到2019年12月31日的死亡率随访数据。
# 公开的数据中仅有随访时间或潜在的死亡原因
# 生存地址原始下载链接：https://ftp.cdc.gov/pub/Health_Statistics/NCHS/datalinkage/linked_mortality/
# 可以下载不同周期的生存数据，下载后进行读取

# setwd("C:/PUBLIC USE DATA")
mort.data_2015_2016 <- NHSMortRead('G:/BaiduNetdiskDownload/NHANES/NHANES_2015_2016_MORT_2019_PUBLIC.dat')
isexit <- file.exists('G:/BaiduNetdiskDownload/NHANES/NHANES_2015_2016_MORT_2019_PUBLIC.dat')
isexit

# eligstat
# A value of 1 for ELIGSTAT indicates that the survey participant was eligible for the mortality linkage, 
# a value of 2 indicates the survey participant was under age 18 and not eligible for public release, 
# a value of 3 indicates the survey participant was not linkage eligible due to having insufficient identifying data to conduct data linkage.
 
# ucod_leading: 10种 潜在死亡原因，001-010，分别对应的原因看课程附件材料
# mortstat: 1生存状态
# permth_int: 从interview到随访时间点，间隔的月份数
# permth_exm: 从MEC-移动检查车到随访时间点，间隔的月份数

mort.data <- mort.data[which(mort.data$eligstat == 1),]
dim(mort.data_2015_2016) # 9971
# 转变为大写的列名，和 demo 进行拼接
colnames(mort.data_2015_2016) <- toupper(colnames(mort.data_2015_2016))
# View(mort.data_2015_2016)


##### 2.2 合并数据提取变量  & 去掉 NA 的行 & 衍生变量 ##### 
analyze.sample.data <- demo.i[,c('SEQN',
                                 "RIDAGEYR", "RIAGENDR", "INDFMPIR", "DMDEDUC2",
                                 "WTINT2YR", "SDMVPSU", "SDMVSTRA")]

analyze.mort.data <- mort.data_2015_2016[,c('SEQN', 'ELIGSTAT', 
                                            "MORTSTAT", "PERMTH_INT")]
# 拼接 demo & mort 数据
analyze.sample.data.add.mort <- plyr::join_all(list(analyze.sample.data, analyze.mort.data),by = 'SEQN')
colnames(analyze.sample.data)
colnames(analyze.mort.data)

dim(analyze.sample.data.add.mort) # 9971    
# View(analyze.sample.data.add.mort) 

# 一键去掉 NA 的行 // 这里一下子删除没了 数据不统一吗? 数据取错了
analyze.sample.data.add.mort.drop.na <- drop_na(analyze.sample.data.add.mort) 
dim(analyze.sample.data.add.mort.drop.na) # 5071

# 衍生低收入 vs 非低收入的变量 PIR.factor, 作为后续演示的 Y 变量
analyze.sample.data.add.mort.drop.na$PIR.factor <- ifelse(analyze.sample.data.add.mort.drop.na$INDFMPIR < 1.3, 1, 0)

##### 2.3 生成复杂抽样 NHANES_design ##### 
NHANES_design <- svydesign(
  data = analyze.sample.data.add.mort.drop.na, 
  ids = ~SDMVPSU, 
  strata = ~SDMVSTRA, 
  nest = TRUE, 
  weights = ~WTINT2YR,
  survey.lonely.psu = "adjust") # 可以加上 survey.lonely.psu = "adjust" 避免1个PSU报错

summary(NHANES_design)

##### 2.4 复杂抽样的 Cox 回归模型 ##### 
# 死亡和贫困指数之间的相关性
m1 <- svycoxph(Surv(PERMTH_INT, MORTSTAT) ~ PIR.factor, 
               design = NHANES_design)
tbl_regression(m1)

m1 <- svycoxph(Surv(PERMTH_INT, MORTSTAT) ~ PIR.factor + RIDAGEYR + factor(DMDEDUC2), 
               design = NHANES_design)

summary(m1)


##### 2.5 使用 gtsummary 优化输出 #####
tbl_regression(m1)

